!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine write_vtu_data(file_unit, var_name, var, x, y, z, nx, ny, nz, n_var, var_dim, format)
!==================================================================================================================================|

   implicit none

   integer(4), parameter :: n_dim = 3
   integer(4), parameter :: n_point_cell = 8

   character(*), intent(in) :: var_name(n_var)
   integer(4), intent(in) :: nx, ny, nz, n_var, format, file_unit

   integer(4), intent(in) :: var_dim(n_var)
   real(4), intent(in) :: x(nx), y(ny), z(nz)
   real(4), intent(in) :: var(sum(var_dim), nx, ny, nz)

   character(12) :: int_num1, int_num2

   integer(4) :: int_size = 4, float_size = 4
   integer(4) :: i_var, i, j, k, index_var
   integer(4), allocatable :: cell_type(:)

   integer(4) :: n_point, n_cell, i_point, pos, i_cell
   integer(4) :: nxp1, nxny, nxnyp1, nxnypnxp1, nxnypnx
   integer(4), allocatable :: connectivity(:, :), offset(:)

   real(4), allocatable :: point(:, :)

!----------------------------------------------------------------------------------------------------------------------------------|
!  prepare
!----------------------------------------------------------------------------------------------------------------------------------|
   n_point = nx * ny * nz
   n_cell = (nx - 1) * (ny - 1) * (nz - 1)
   nxp1 = nx + 1
   nxny = nx * ny
   nxnyp1 = nx * ny + 1
   nxnypnxp1 = nx * ny + nx + 1
   nxnypnx = nx * ny + nx

   allocate(cell_type(n_cell))
   allocate(connectivity(n_point_cell, n_cell))
   allocate(offset(n_cell))
   allocate(point(n_dim, n_point))
   cell_type = 12

   i_point = 1
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      point(1, i_point) = x(i)
      point(2, i_point) = y(j)
      point(3, i_point) = z(k)
      i_point = i_point + 1
   enddo
   enddo
   enddo

   i_point = 1
   do k = 1, nz - 1
   do j = 1, ny - 1
   do i = 1, nx - 1
      pos = (i - 1) * 1 + (j - 1) * nx + (k - 1) * nxny
      connectivity(1, i_point) = pos
      connectivity(2, i_point) = pos + 1
      connectivity(3, i_point) = pos + nxp1
      connectivity(4, i_point) = pos + nx
      connectivity(5, i_point) = pos + nxny
      connectivity(6, i_point) = pos + nxnyp1
      connectivity(7, i_point) = pos + nxnypnxp1
      connectivity(8, i_point) = pos + nxnypnx
      i_point = i_point + 1
   enddo
   enddo
   enddo

   do i_cell = 1, n_cell
      offset(i_cell) = i_cell * n_point_cell
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  binary
!----------------------------------------------------------------------------------------------------------------------------------|
   if (format .eq. 0) then
      index_var = 1
      do i_var = 1, n_var
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the scalar point data
!----------------------------------------------------------------------------------------------------------------------------------|
         if (var_dim(i_var) .eq. 1) then
            write(file_unit) n_point * float_size, var(index_var, :, :, :)
            index_var = index_var + var_dim(i_var)
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the vector point data
!----------------------------------------------------------------------------------------------------------------------------------|
         elseif (var_dim(i_var) .eq. 3) then
            write(file_unit) 3 * n_point * float_size, var(index_var : index_var + 2, :, :, :)
            index_var = index_var + var_dim(i_var)
         endif
      enddo
      write(file_unit) 3 * n_point * float_size, point
      write(file_unit) n_point_cell * n_cell * int_size, connectivity
      write(file_unit) n_cell * int_size, offset
      write(file_unit) n_cell * int_size, cell_type

!----------------------------------------------------------------------------------------------------------------------------------|
!  ASCII
!----------------------------------------------------------------------------------------------------------------------------------|
   else
!----------------------------------------------------------------------------------------------------------------------------------|
!     write the info data
!----------------------------------------------------------------------------------------------------------------------------------|
      write(int_num1, '(i12)') n_point
      write(int_num2, '(i12)') n_cell
      write(file_unit, '(a)') '<Piece NumberOfPoints="' // int_num1 // '" NumberOfCells="' // int_num2 // '">'

      write(file_unit, '(a)') '<PointData>'
      index_var = 1
      do i_var = 1, n_var
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the scalar point data
!----------------------------------------------------------------------------------------------------------------------------------|
         if (var_dim(i_var) .eq. 1) then
            write(file_unit, '(a)') '<DataArray type="Float32" format="ascii" NumberOfComponents="1" Name="' //                    &
               trim(var_name(i_var)) // '">'
            write(file_unit, '(6es)') var(index_var, :, :, :)
            write(file_unit, '(a)') '</DataArray>'
            index_var = index_var + var_dim(i_var)
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the vector point data
!----------------------------------------------------------------------------------------------------------------------------------|
         elseif (var_dim(i_var) .eq. 3) then
            write(file_unit, '(a)') '<DataArray type="Float32" format="ascii" NumberOfComponents="3" Name="' //                 &
               trim(var_name(i_var)) // '">'
            write(file_unit, '(6es)') var(index_var : index_var + 2, :, :, :)
            write(file_unit, '(a)') '</DataArray>'
            index_var = index_var + var_dim(i_var)
         endif
      enddo
      write(file_unit, '(a)') '</PointData>'

      write(file_unit, '(a)') '<Points>'
      write(file_unit, '(a)') '<DataArray type="Float32" format="ascii" Name="coordinates" NumberOfComponents="3">'
      write(file_unit, '(6es)') point
      write(file_unit, '(a)') '</DataArray>'
      write(file_unit, '(a)') '</Points>'

      write(file_unit, '(a)') '<Cells>'
      write(file_unit, '(a)') '<DataArray type="Int32" format="ascii" Name="connectivity">'
      write(file_unit, '(9i)') connectivity
      write(file_unit, '(a)') '</DataArray>'
      write(file_unit, '(a)') '<DataArray type="Int32" format="ascii" Name="offsets">'
      write(file_unit, '(12i)') offset
      write(file_unit, '(a)') '</DataArray>'
      write(file_unit, '(a)') '<DataArray type="Int32" format="ascii" Name="types">'
      write(file_unit, '(12i)') cell_type
      write(file_unit, '(a)') '</DataArray>'
      write(file_unit, '(a)') '</Cells>'
      write(file_unit, '(a)') '</Piece>'
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  free memory
!----------------------------------------------------------------------------------------------------------------------------------|
   deallocate (cell_type, connectivity, offset, point)

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine write_vtu_data
